1 Dapprima selezione delle variabili per ciascuna delle quattro specie

2 Data

Thanks Marco Salvatori.

covar<- read.csv2('data/covs_original_ALL.csv', stringsAsFactors = F)
dt1<- covar[,c('Sampling.Unit', 'Area', 'CT.sens', 'EL', 'SL', 'Dis')]
dt1[,4:6] <- apply(X=dt1[,4:6], MARGIN=2, FUN=decostand, method='standardize')
var<- dt1
ll<-readRDS(file = 'data/matrici_complessive.rds')
dom1<-ll$domestici
dom<-ll$domestici
ibex<-ll$ibex
leo<-ll$leo
lupus<-ll$lupo

3 Single species

DOMESTICI

do<- unmarkedFrameOccu(y=dom, siteCovs = var, obsCovs = NULL)
ib<- unmarkedFrameOccu(y=ibex, siteCovs = var, obsCovs = NULL)
le<- unmarkedFrameOccu(y=leo, siteCovs = var, obsCovs = NULL)
lu<- unmarkedFrameOccu(y=lupus, siteCovs = var, obsCovs = NULL)

3.1 Bestiame

##p
mod0<-occu(~1~Area+EL+SL+Dis ,do)
mod1<-occu(~Dis~Area+EL+SL+Dis,do)
mod2<-occu(~CT.sens~Area+EL+SL+Dis,do)
mod3<-occu(~Dis+CT.sens~Area+EL+SL+Dis, do)
mods<-fitList('Zero'=mod0,'Dis'=mod1,'CT.sens'=mod2, 'Dis+CT.sens'=mod3)
modSel(mods)
## Hessian is singular.
##             nPars     AIC  delta    AICwt cumltvWt
## Dis+CT.sens    10 5230.07   0.00  1.0e+00     1.00
## CT.sens         9 5317.00  86.93  1.3e-19     1.00
## Zero            8 5363.11 133.03  1.3e-29     1.00
## Dis             9 5832.17 602.09 1.8e-131     1.00
## psi
mod1<-occu(~Dis+CT.sens~1, do)
mod2<-occu(~Dis+CT.sens~Area, do)
mod3<-occu(~Dis+CT.sens~Dis, do)
mod4<-occu(~Dis+CT.sens~EL, do)
mod5<-occu(~Dis+CT.sens~SL, do)
mod6<-occu(~Dis+CT.sens~Area+SL, do)
mod7<-occu(~Dis+CT.sens~Area+EL, do)
mod8<-occu(~Dis+CT.sens~Area+Dis, do)
mod9<-occu(~Dis+CT.sens~EL+Dis, do)
mod10<-occu(~Dis+CT.sens~EL+SL, do)
mod11<-occu(~Dis+CT.sens~Dis+SL, do)
mod12<-occu(~Dis+CT.sens~Area+EL+SL, do)
mod13<-occu(~Dis+CT.sens~Area+EL+Dis, do)
mod14<-occu(~Dis+CT.sens~Area+SL+Dis, do)
mod15<-occu(~Dis+CT.sens~SL+EL+Dis, do)
mod16<-occu(~Dis+CT.sens~ Area+EL+SL+Dis, do)
mdls<-fitList('Zero'=mod0,'K'=mod1,'Area'=mod2,'Dis'=mod3,
              'EL'=mod4,'SL'=mod5,'Area+SL'=mod6,'Area+EL'=mod7, 
              'Area+Dis'=mod8, 'EL+Dis'=mod9, 'EL+SL'=mod10, 'Dis+SL'=mod11,
              'Area+EL+SL'=mod12, 'Area+EL+Dis'=mod13,  'Area+SL+Dis'=mod14, 
              'SL+EL+Dis'=mod15, 'Area+EL+SL+Dis'=mod16)
ms<-modSel(mdls)
msms<- ms@Full[,c('model', 'negLogLike', 'nPars', 'AIC', 'delta', 'AICwt')]
msms
##             model negLogLike nPars      AIC      delta         AICwt
## 14    Area+EL+Dis   2605.266     9 5228.533   0.000000  6.755202e-01
## 17 Area+EL+SL+Dis   2605.036    10 5230.072   1.539166  3.129051e-01
## 8         Area+EL   2611.004     8 5238.007   9.474579  5.919155e-03
## 13     Area+EL+SL   2610.830     9 5239.660  11.127025  2.590810e-03
## 10         EL+Dis   2614.273     6 5240.545  12.012452  1.664055e-03
## 16      SL+EL+Dis   2614.034     7 5242.069  13.536115  7.767984e-04
## 5              EL   2616.613     5 5243.226  14.693349  4.355301e-04
## 11          EL+SL   2616.469     6 5244.939  16.406348  1.849461e-04
## 3            Area   2619.938     7 5253.875  25.342434  2.121283e-06
## 7         Area+SL   2619.828     8 5255.657  27.123998  8.704361e-07
## 12         Dis+SL   2622.515     6 5257.029  28.496660  4.381948e-07
## 2               K   2627.666     4 5263.332  34.799490  1.875103e-08
## 6              SL   2627.600     5 5265.199  36.666602  7.372025e-09
## 1            Zero   2673.553     8 5363.106 134.572937  4.050386e-30
## 15    Area+SL+Dis   2868.362     9 5754.725 526.192226 3.702100e-115
## 4             Dis   2904.134     5 5818.269 589.736298 5.888821e-129
## 9        Area+Dis   2904.189     8 5824.379 595.846176 2.775142e-130
summary(mod13)
## 
## Call:
## occu(formula = ~Dis + CT.sens ~ Area + EL + Dis, data = do)
## 
## Occupancy (logit-scale):
##             Estimate    SE      z P(>|z|)
## (Intercept)    0.131 0.277  0.471 0.63789
## AreaSB        -0.688 0.416 -1.654 0.09803
## AreaSU         1.008 0.475  2.121 0.03389
## AreaTB        -0.811 0.459 -1.765 0.07758
## EL            -0.582 0.187 -3.111 0.00187
## Dis           -0.617 0.190 -3.241 0.00119
## 
## Detection (logit-scale):
##               Estimate     SE      z  P(>|z|)
## (Intercept)     -2.436 0.0504 -48.34 0.00e+00
## Dis             -0.467 0.0552  -8.47 2.53e-17
## CT.sensmedium    0.487 0.0850   5.73 1.02e-08
## 
## AIC: 5228.533 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 87 
## Bootstrap iterations: 0

3.2 Stambecco Siberiano

### p

mod0<-occu(~1~Area+EL+SL+Dis ,ib)
mod1<-occu(~Dis~Area+EL+SL+Dis,ib)
mod2<-occu(~CT.sens~Area+EL+SL+Dis,ib)
mod3<-occu(~Dis+CT.sens~Area+EL+SL+Dis, ib)
mods<-fitList('Zero'=mod0,'Dis'=mod1,'CT.sens'=mod2, 'Dis+CT.sens'=mod3)
modSel(mods)
##             nPars     AIC  delta   AICwt cumltvWt
## Dis+CT.sens    10 1444.34   0.00 9.6e-01     0.96
## Dis             9 1450.47   6.13 4.5e-02     1.00
## Zero            8 1558.03 113.68 2.0e-25     1.00
## CT.sens         9 1560.50 116.16 5.7e-26     1.00
## psi
mod1<-occu(~CT.sens~1, ib)
mod2<-occu(~CT.sens~Area, ib)
mod3<-occu(~CT.sens~Dis, ib)
mod4<-occu(~CT.sens~EL, ib)
mod5<-occu(~CT.sens~SL, ib)
mod6<-occu(~CT.sens~Area+SL, ib)
mod7<-occu(~CT.sens~Area+EL, ib)
mod8<-occu(~CT.sens~Area+Dis, ib)
mod9<-occu(~CT.sens~EL+Dis, ib)
mod10<-occu(~CT.sens~EL+SL, ib)
mod11<-occu(~CT.sens~Dis+SL, ib)
mod12<-occu(~CT.sens~Area+EL+SL, ib)
mod13<-occu(~CT.sens~Area+EL+Dis, ib)
mod14<-occu(~CT.sens~Area+SL+Dis, ib)
mod15<-occu(~CT.sens~SL+EL+Dis, ib)
mod16<-occu(~CT.sens~ Area+EL+SL+Dis, ib)
mdls<-fitList('Zero'=mod0,'K'=mod1,'Area'=mod2,'Dis'=mod3,
              'EL'=mod4,'SL'=mod5,'Area+SL'=mod6,'Area+EL'=mod7, 
              'Area+Dis'=mod8, 'EL+Dis'=mod9, 'EL+SL'=mod10, 'Dis+SL'=mod11,
              'Area+EL+SL'=mod12, 'Area+EL+Dis'=mod13,  'Area+SL+Dis'=mod14, 
              'SL+EL+Dis'=mod15, 'Area+EL+SL+Dis'=mod16)
ms<-modSel(mdls)
## Hessian is singular.
## Hessian is singular.
msms<- ms@Full[,c('model', 'negLogLike', 'nPars', 'AIC', 'delta', 'AICwt')]
msms
##             model negLogLike nPars      AIC      delta        AICwt
## 13     Area+EL+SL   714.4724     8 1444.945   0.000000 4.711932e-01
## 6              SL   719.1642     4 1446.328   1.383615 2.359124e-01
## 11          EL+SL   718.6172     5 1447.234   2.289506 1.499820e-01
## 16      SL+EL+Dis   718.3004     6 1448.601   3.656000 7.573712e-02
## 5              EL   721.3001     4 1450.600   5.655432 2.787004e-02
## 4             Dis   721.3584     4 1450.717   5.771906 2.629333e-02
## 10         EL+Dis   721.0618     5 1452.124   7.178808 1.301190e-02
## 8         Area+EL   751.9379     7 1517.876  72.930922 6.861998e-17
## 3            Area   753.7809     6 1519.562  74.617022 2.953375e-17
## 14    Area+EL+Dis   751.9234     8 1519.847  74.901928 2.561251e-17
## 9        Area+Dis   753.6673     7 1521.335  76.389846 1.217180e-17
## 1            Zero   771.0131     8 1558.026 113.081478 1.311801e-25
## 17 Area+EL+SL+Dis   771.2512     9 1560.502 115.557665 3.803390e-26
## 15    Area+SL+Dis   778.7924     8 1573.585 128.640099 5.487271e-29
## 7         Area+SL   781.5926     7 1577.185 132.240404 9.069018e-30
## 2               K   794.2281     3 1594.456 149.511353 1.611467e-33
## 12         Dis+SL   794.2286     5 1598.457 153.512482 2.179652e-34
summary(mod12)
## 
## Call:
## occu(formula = ~CT.sens ~ Area + EL + SL, data = ib)
## 
## Occupancy (logit-scale):
##             Estimate    SE      z P(>|z|)
## (Intercept)   -0.682 0.307 -2.221  0.0263
## AreaSB        -0.576 0.490 -1.175  0.2399
## AreaSU         0.174 0.424  0.411  0.6812
## AreaTB        -1.295 0.615 -2.106  0.0352
## EL             0.270 0.187  1.448  0.1476
## SL             0.368 0.169  2.173  0.0298
## 
## Detection (logit-scale):
##               Estimate    SE      z   P(>|z|)
## (Intercept)     -3.520 0.112 -31.55 1.92e-218
## CT.sensmedium    0.684 0.227   3.02  2.55e-03
## 
## AIC: 1444.945 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 130 
## Bootstrap iterations: 0

3.3 Leopardo delle nevi

### p

mod0<-occu(~1~Area+EL+SL+Dis ,le)
mod1<-occu(~Dis~Area+EL+SL+Dis,le)
mod2<-occu(~CT.sens~Area+EL+SL+Dis,le)
mod3<-occu(~Dis+CT.sens~Area+EL+SL+Dis, le)
mods<-fitList('Zero'=mod0,'Dis'=mod1,'CT.sens'=mod2, 'Dis+CT.sens'=mod3)#, 'Dis+CT.sens+Area'=mod4 )
modSel(mods)
##             nPars     AIC delta   AICwt cumltvWt
## Dis+CT.sens    10 2481.82  0.00 5.5e-01     0.55
## CT.sens         9 2482.22  0.41 4.5e-01     1.00
## Dis             9 2524.36 42.54 3.2e-10     1.00
## Zero            8 2531.31 49.50 9.8e-12     1.00
## psi
mod1<-occu(~CT.sens~1, le)
mod2<-occu(~CT.sens~Area, le, starts = rep(-1,6))
mod3<-occu(~CT.sens~Dis, le)
mod4<-occu(~CT.sens~EL, le)
mod5<-occu(~CT.sens~SL, le)
mod6<-occu(~CT.sens~Area+SL, le, starts = rep(-1,7))
mod7<-occu(~CT.sens~Area+EL, le)
mod8<-occu(~CT.sens~Area+Dis, le, starts = rep(1,7))
mod9<-occu(~CT.sens~EL+Dis, le)
mod10<-occu(~CT.sens~EL+SL, le)
mod11<-occu(~CT.sens~Dis+SL, le)
mod12<-occu(~CT.sens~Area+EL+SL, le)
mod13<-occu(~CT.sens~Area+EL+Dis, le, starts = rep(1,8))
mod14<-occu(~CT.sens~Area+SL+Dis, le, starts = rep(1,8))
mod15<-occu(~CT.sens~SL+EL+Dis, le)
mod16<-occu(~CT.sens~ Area+EL+SL+Dis, le)
mdls<-fitList('Zero'=mod0,'K'=mod1,'Area'=mod2,'Dis'=mod3,
              'EL'=mod4,'SL'=mod5,'Area+SL'=mod6,'Area+EL'=mod7, 
              'Area+Dis'=mod8, 'EL+Dis'=mod9, 'EL+SL'=mod10, 'Dis+SL'=mod11,
              'Area+EL+SL'=mod12, 'Area+EL+Dis'=mod13,  'Area+SL+Dis'=mod14, 
              'SL+EL+Dis'=mod15, 'Area+EL+SL+Dis'=mod16)
ms<-modSel(mdls)
msms<- ms@Full[,c('model', 'negLogLike', 'nPars', 'AIC', 'delta', 'AICwt')]
msms
##             model negLogLike nPars      AIC       delta        AICwt
## 3            Area   1233.772     6 2479.543  0.00000000 2.089589e-01
## 9        Area+Dis   1232.798     7 2479.596  0.05234249 2.035611e-01
## 7         Area+SL   1232.951     7 2479.903  0.35945726 1.745845e-01
## 15    Area+SL+Dis   1232.205     8 2480.410  0.86671176 1.354742e-01
## 14    Area+EL+Dis   1232.715     8 2481.429  1.88585873 8.138639e-02
## 8         Area+EL   1233.769     7 2481.539  1.99546759 7.704607e-02
## 13     Area+EL+SL   1232.951     8 2481.903  2.35929613 6.423122e-02
## 17 Area+EL+SL+Dis   1232.111     9 2482.222  2.67844085 5.475764e-02
## 10         EL+Dis   1257.842     5 2525.685 46.14126918 1.998074e-11
## 16      SL+EL+Dis   1257.513     6 2527.026 47.48274234 1.021679e-11
## 1            Zero   1257.656     8 2531.312 51.76872549 1.198464e-12
## 4             Dis   1262.191     4 2532.382 52.83916398 7.017495e-13
## 12         Dis+SL   1261.839     5 2533.678 54.13488590 3.671298e-13
## 5              EL   1267.189     4 2542.378 62.83454751 4.739278e-15
## 2               K   1268.284     3 2542.569 63.02572991 4.307224e-15
## 11          EL+SL   1266.984     5 2543.968 64.42488156 2.139811e-15
## 6              SL   1268.066     4 2544.132 64.58845437 1.971769e-15
summary(mod2)
## 
## Call:
## occu(formula = ~CT.sens ~ Area, data = le, starts = rep(-1, 6))
## 
## Occupancy (logit-scale):
##             Estimate    SE     z  P(>|z|)
## (Intercept)   -0.753 0.305 -2.47 1.36e-02
## AreaSB         1.805 1.599  1.13 2.59e-01
## AreaSU         2.293 0.465  4.93 8.36e-07
## AreaTB        -2.139 0.789 -2.71 6.69e-03
## 
## Detection (logit-scale):
##               Estimate     SE      z  P(>|z|)
## (Intercept)      -3.05 0.0661 -46.16 0.000000
## CT.sensmedium    -1.57 0.4430  -3.54 0.000405
## 
## AIC: 2479.543 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 26 
## Bootstrap iterations: 0

3.4 Lupo

### p

mod0<-occu(~1~Area+EL+SL+Dis ,lu)
mod1<-occu(~Dis~Area+EL+SL+Dis,lu)
mod2<-occu(~CT.sens~Area+EL+SL+Dis,lu)
mod3<-occu(~Dis+CT.sens~Area+EL+SL+Dis, lu)
mods<-fitList('Zero'=mod0,'Dis'=mod1,'CT.sens'=mod2, 'Dis+CT.sens'=mod3)#, 'Dis+CT.sens+Area'=mod4 )
modSel(mods)
##             nPars     AIC delta AICwt cumltvWt
## Zero            8 1043.98  0.00 0.523     0.52
## CT.sens         9 1045.83  1.85 0.208     0.73
## Dis             9 1045.98  2.00 0.193     0.92
## Dis+CT.sens    10 1047.82  3.84 0.077     1.00
## psi
mod1<-occu(~1~1, lu)
mod2<-occu(~1~Area, lu)
mod3<-occu(~1~Dis, lu)
mod4<-occu(~1~EL, lu)
mod5<-occu(~1~SL, lu)
mod6<-occu(~1~Area+SL, lu)
mod7<-occu(~1~Area+EL, lu)
mod8<-occu(~1~Area+Dis, lu)
mod9<-occu(~1~EL+Dis, lu)
mod10<-occu(~1~EL+SL, lu)
mod11<-occu(~1~Dis+SL, lu)
mod12<-occu(~1~Area+EL+SL, lu)
mod13<-occu(~1~Area+EL+Dis, lu)
mod14<-occu(~1~Area+SL+Dis, lu)
mod15<-occu(~1~SL+EL+Dis, lu)
mod16<-occu(~1~ Area+EL+SL+Dis, lu)
mdls<-fitList('Zero'=mod0,'K'=mod1,'Area'=mod2,'Dis'=mod3,
              'EL'=mod4,'SL'=mod5,'Area+SL'=mod6,'Area+EL'=mod7, 
              'Area+Dis'=mod8, 'EL+Dis'=mod9, 'EL+SL'=mod10, 'Dis+SL'=mod11,
              'Area+EL+SL'=mod12, 'Area+EL+Dis'=mod13,  'Area+SL+Dis'=mod14, 
              'SL+EL+Dis'=mod15, 'Area+EL+SL+Dis'=mod16)
ms<-modSel(mdls)
msms<- ms@Full[,c('model', 'negLogLike', 'nPars', 'AIC', 'delta', 'AICwt')]
msms
##             model negLogLike nPars      AIC    delta       AICwt
## 10         EL+Dis   515.8354     4 1039.671 0.000000 0.254151574
## 4             Dis   516.9341     3 1039.868 0.197321 0.230274105
## 16      SL+EL+Dis   515.7615     5 1041.523 1.852063 0.100675192
## 12         Dis+SL   516.8750     4 1041.750 2.079139 0.089869740
## 14    Area+EL+Dis   514.1072     7 1042.214 2.543562 0.071246811
## 2               K   519.6448     2 1043.290 3.618671 0.041620605
## 9        Area+Dis   515.7373     6 1043.475 3.803640 0.037943976
## 1            Zero   513.9901     8 1043.980 4.309232 0.029468296
## 17 Area+EL+SL+Dis   513.9901     8 1043.980 4.309232 0.029468296
## 3            Area   517.1861     5 1044.372 4.701246 0.024223126
## 5              EL   519.3935     3 1044.787 5.116011 0.019686343
## 8         Area+EL   516.6010     6 1045.202 5.531086 0.015996782
## 15    Area+SL+Dis   515.6164     7 1045.233 5.561959 0.015751745
## 6              SL   519.6207     3 1045.241 5.570524 0.015684436
## 7         Area+SL   517.0778     6 1046.156 6.484662 0.009930393
## 11          EL+SL   519.3674     4 1046.735 7.064001 0.007433007
## 13     Area+EL+SL   516.4900     7 1046.980 7.309139 0.006575572
summary(mod9)
## 
## Call:
## occu(formula = ~1 ~ EL + Dis, data = lu)
## 
## Occupancy (logit-scale):
##             Estimate    SE     z  P(>|z|)
## (Intercept)   -0.764 0.221 -3.46 0.000535
## EL            -0.299 0.205 -1.46 0.143235
## Dis            0.506 0.198  2.55 0.010718
## 
## Detection (logit-scale):
##  Estimate    SE     z   P(>|z|)
##     -3.95 0.143 -27.6 8.16e-168
## 
## AIC: 1039.671 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 49 
## Bootstrap iterations: 0

3.5 Synthesis

Perciò i migliori modelli di singola specie, selezionando il migliore in termini di AIC, sono i seguenti:

##                          p                      psi
## leopardo          ~CT.sens                    ~Area
## lupo                    ~1      ~Elevation+Distance
## stambecco         ~CT.sens    ~Area+Elevation+Slope
## bestiame  Distance+CT.sens ~Area+Elevation+Distance

3.6 Occupancy di singola specie - grafici

4 Test delle ipotesi ecologiche di interazione inter-specifica

4.1 Interazione stambecco - bestiame domestico

dom_ibex_dip <- occuMulti(detformulas = c('~Dis+CT.sens','~CT.sens'), stateformulas = c('~Area+EL+Dis', '~Area+EL+SL', '~1'), data = dt1, control=list(maxit=1000, trace=TRUE, REPORT=1))
dom_ibex_ind <- occuMulti(detformulas = c('~Dis+CT.sens','~CT.sens'), stateformulas = c('~Area+EL+Dis', '~Area+EL+SL', '~0'), data = dt1, control=list(maxit=1000, trace=TRUE, REPORT=1))
library(AICcmodavg)
modavgShrink(list(dom_ibex_dip, dom_ibex_ind), parm='[domestici:ibex] (Intercept)', parm.type='psi')
##          model negLogLike nPars      AIC      delta     AICwt
## 2 dom_ibex_ind   3319.739    17 6673.477 0.00000000 0.5093794
## 1 dom_ibex_dip   3318.776    18 6673.552 0.07504392 0.4906206

Piccola differenza, ma il modello con l’AIC migliore include l’interazione domestici-stambecco.

Modello senza interazione:

## 
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~CT.sens"), stateformulas = c("~Area+EL+Dis", 
##     "~Area+EL+SL", "~0"), data = dt1, control = list(maxit = 1000, 
##     trace = TRUE, REPORT = 1), maxOrder = 2L)
## 
## Occupancy (logit-scale):
##                         Estimate    SE      z P(>|z|)
## [domestici] (Intercept)    0.130 0.277  0.471 0.63796
## [domestici] AreaSB        -0.687 0.416 -1.654 0.09816
## [domestici] AreaSU         1.008 0.475  2.121 0.03389
## [domestici] AreaTB        -0.811 0.459 -1.766 0.07748
## [domestici] EL            -0.582 0.187 -3.111 0.00187
## [domestici] Dis           -0.617 0.190 -3.241 0.00119
## [ibex] (Intercept)        -0.681 0.307 -2.220 0.02639
## [ibex] AreaSB             -0.576 0.490 -1.175 0.23987
## [ibex] AreaSU              0.174 0.424  0.411 0.68141
## [ibex] AreaTB             -1.295 0.615 -2.106 0.03517
## [ibex] EL                  0.270 0.187  1.448 0.14756
## [ibex] SL                  0.368 0.169  2.173 0.02980
## 
## Detection (logit-scale):
##                           Estimate     SE      z   P(>|z|)
## [domestici] (Intercept)     -2.436 0.0504 -48.34  0.00e+00
## [domestici] Dis             -0.467 0.0552  -8.47  2.53e-17
## [domestici] CT.sensmedium    0.487 0.0850   5.73  1.02e-08
## [ibex] (Intercept)          -3.520 0.1116 -31.55 1.97e-218
## [ibex] CT.sensmedium         0.684 0.2266   3.02  2.54e-03
## 
## AIC: 6673.477 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 88 
## Bootstrap iterations: 0

Modello con interazione:

## 
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~CT.sens"), stateformulas = c("~Area+EL+Dis", 
##     "~Area+EL+SL", "~1"), data = dt1, control = list(maxit = 1000, 
##     trace = TRUE, REPORT = 1), maxOrder = 2L)
## 
## Occupancy (logit-scale):
##                              Estimate    SE      z P(>|z|)
## [domestici] (Intercept)         0.318 0.312  1.020 0.30796
## [domestici] AreaSB             -0.766 0.423 -1.809 0.07049
## [domestici] AreaSU              1.027 0.479  2.146 0.03190
## [domestici] AreaTB             -0.938 0.472 -1.988 0.04682
## [domestici] EL                 -0.564 0.189 -2.992 0.00277
## [domestici] Dis                -0.616 0.190 -3.239 0.00120
## [ibex] (Intercept)             -0.398 0.370 -1.076 0.28198
## [ibex] AreaSB                  -0.687 0.501 -1.372 0.16994
## [ibex] AreaSU                   0.234 0.430  0.544 0.58620
## [ibex] AreaTB                  -1.377 0.620 -2.221 0.02632
## [ibex] EL                       0.184 0.198  0.930 0.35230
## [ibex] SL                       0.366 0.169  2.169 0.03006
## [domestici:ibex] (Intercept)   -0.522 0.379 -1.376 0.16867
## 
## Detection (logit-scale):
##                           Estimate     SE      z   P(>|z|)
## [domestici] (Intercept)     -2.437 0.0504 -48.32  0.00e+00
## [domestici] Dis             -0.467 0.0552  -8.47  2.53e-17
## [domestici] CT.sensmedium    0.487 0.0850   5.73  1.01e-08
## [ibex] (Intercept)          -3.522 0.1118 -31.51 6.95e-218
## [ibex] CT.sensmedium         0.687 0.2263   3.04  2.39e-03
## 
## AIC: 6673.552 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 152 
## Bootstrap iterations: 0

4.2 Interazione leopardo - bestiame domestico

##         model negLogLike nPars      AIC    delta     AICwt
## 1 dom_leo_dip   3837.190    16 7706.380 0.000000 0.7496418
## 2 dom_leo_ind   3839.286    15 7708.573 2.193406 0.2503582

Il modello migliore per AIC contiene l’interazione domestici-leopardo.

Modello senza interazione:

## 
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~CT.sens"), stateformulas = c("~Area+EL+Dis", 
##     "~Area", "~0"), data = dt2, control = list(maxit = 1000, 
##     trace = TRUE, REPORT = 1), maxOrder = 2L)
## 
## Occupancy (logit-scale):
##                         Estimate    SE      z  P(>|z|)
## [domestici] (Intercept)    0.130 0.277  0.469 6.39e-01
## [domestici] AreaSB        -0.689 0.416 -1.657 9.74e-02
## [domestici] AreaSU         1.009 0.475  2.124 3.37e-02
## [domestici] AreaTB        -0.810 0.459 -1.763 7.79e-02
## [domestici] EL            -0.582 0.187 -3.112 1.86e-03
## [domestici] Dis           -0.616 0.190 -3.239 1.20e-03
## [leopardo] (Intercept)    -0.724 0.303 -2.388 1.70e-02
## [leopardo] AreaSB          9.017   NaN    NaN      NaN
## [leopardo] AreaSU          2.263 0.464  4.877 1.08e-06
## [leopardo] AreaTB         -2.167 0.789 -2.746 6.04e-03
## 
## Detection (logit-scale):
##                           Estimate     SE      z  P(>|z|)
## [domestici] (Intercept)     -2.436 0.0504 -48.34 0.00e+00
## [domestici] Dis             -0.467 0.0552  -8.47 2.51e-17
## [domestici] CT.sensmedium    0.487 0.0850   5.73 1.03e-08
## [leopardo] (Intercept)      -3.052 0.0661 -46.15 0.00e+00
## [leopardo] CT.sensmedium    -1.840 0.2427  -7.58 3.39e-14
## 
## AIC: 7708.573 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 75 
## Bootstrap iterations: 0

Modello con interazione:

## 
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~CT.sens"), stateformulas = c("~Area+EL+Dis", 
##     "~Area", "~1"), data = dt2, control = list(maxit = 1000, 
##     trace = TRUE, REPORT = 1), maxOrder = 2L)
## 
## Occupancy (logit-scale):
##                                  Estimate     SE      z  P(>|z|)
## [domestici] (Intercept)            0.4706  0.340  1.383 1.67e-01
## [domestici] AreaSB                -0.0494  0.588 -0.084 9.33e-01
## [domestici] AreaSU                 1.4974  0.548  2.731 6.31e-03
## [domestici] AreaTB                -1.1095  0.495 -2.242 2.50e-02
## [domestici] EL                    -0.5819  0.187 -3.110 1.87e-03
## [domestici] Dis                   -0.6140  0.190 -3.232 1.23e-03
## [leopardo] (Intercept)            -0.1707  0.418 -0.408 6.83e-01
## [leopardo] AreaSB                  6.2365 54.641  0.114 9.09e-01
## [leopardo] AreaSU                  2.4152  0.491  4.919 8.70e-07
## [leopardo] AreaTB                 -2.5134  0.821 -3.061 2.21e-03
## [domestici:leopardo] (Intercept)  -0.9819  0.507 -1.935 5.30e-02
## 
## Detection (logit-scale):
##                           Estimate     SE      z  P(>|z|)
## [domestici] (Intercept)     -2.437 0.0505 -48.30 0.00e+00
## [domestici] Dis             -0.467 0.0552  -8.47 2.51e-17
## [domestici] CT.sensmedium    0.487 0.0850   5.73 1.02e-08
## [leopardo] (Intercept)      -3.053 0.0661 -46.16 0.00e+00
## [leopardo] CT.sensmedium    -1.834 0.3154  -5.81 6.06e-09
## 
## AIC: 7706.38 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 85 
## Bootstrap iterations: 0

4.3 Interazione lupo - bestiame domestico

##         model negLogLike nPars      AIC    delta     AICwt
## 1 dom_lup_dip   3118.852    14 6265.704 0.000000 0.7772587
## 2 dom_lup_ind   3121.102    13 6268.203 2.499525 0.2227413

Anche in questo caso il modello con migliore AIC prevede l’interazione.

Modello senza interazione:

## 
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~1"), stateformulas = c("~Area+EL+Dis", 
##     "~EL+Dis", "~0"), data = dt3, control = list(maxit = 1000, 
##     trace = TRUE, REPORT = 1), maxOrder = 2L)
## 
## Occupancy (logit-scale):
##                         Estimate    SE      z  P(>|z|)
## [domestici] (Intercept)    0.131 0.277  0.471 0.637865
## [domestici] AreaSB        -0.688 0.416 -1.655 0.098020
## [domestici] AreaSU         1.008 0.475  2.121 0.033895
## [domestici] AreaTB        -0.811 0.459 -1.765 0.077574
## [domestici] EL            -0.582 0.187 -3.111 0.001867
## [domestici] Dis           -0.617 0.190 -3.241 0.001190
## [lupo] (Intercept)        -0.764 0.221 -3.462 0.000537
## [lupo] EL                 -0.299 0.205 -1.463 0.143350
## [lupo] Dis                 0.506 0.198  2.552 0.010715
## 
## Detection (logit-scale):
##                           Estimate     SE      z   P(>|z|)
## [domestici] (Intercept)     -2.436 0.0504 -48.34  0.00e+00
## [domestici] Dis             -0.467 0.0552  -8.47  2.53e-17
## [domestici] CT.sensmedium    0.487 0.0850   5.73  1.02e-08
## [lupo] (Intercept)          -3.950 0.1431 -27.61 8.46e-168
## 
## AIC: 6268.203 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 118 
## Bootstrap iterations: 0

Modello con interazione:

## 
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~1"), stateformulas = c("~Area+EL+Dis", 
##     "~EL+Dis", "~1"), data = dt3, control = list(maxit = 1000, 
##     trace = TRUE, REPORT = 1), maxOrder = 2L)
## 
## Occupancy (logit-scale):
##                              Estimate    SE      z  P(>|z|)
## [domestici] (Intercept)        -0.154 0.310 -0.495 0.620390
## [domestici] AreaSB             -0.688 0.416 -1.656 0.097648
## [domestici] AreaSU              1.002 0.475  2.111 0.034794
## [domestici] AreaTB             -0.813 0.460 -1.768 0.077043
## [domestici] EL                 -0.548 0.191 -2.873 0.004071
## [domestici] Dis                -0.733 0.205 -3.579 0.000345
## [lupo] (Intercept)             -1.300 0.341 -3.807 0.000140
## [lupo] EL                      -0.187 0.209 -0.893 0.372037
## [lupo] Dis                      0.599 0.210  2.848 0.004406
## [domestici:lupo] (Intercept)    0.912 0.441  2.068 0.038680
## 
## Detection (logit-scale):
##                           Estimate     SE      z   P(>|z|)
## [domestici] (Intercept)     -2.436 0.0504 -48.35  0.00e+00
## [domestici] Dis             -0.467 0.0552  -8.47  2.55e-17
## [domestici] CT.sensmedium    0.487 0.0850   5.72  1.05e-08
## [lupo] (Intercept)          -3.939 0.1418 -27.77 1.06e-169
## 
## AIC: 6265.704 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 144 
## Bootstrap iterations: 0

4.4 Interazione leopardo - stambecco

##        model negLogLike nPars      AIC    delta     AICwt
## 1 ib_leo_dip   1944.109    15 3918.218 0.000000 0.9583079
## 2 ib_leo_ind   1948.244    14 3924.488 6.269715 0.0416921

In questo caso il modello con AIC più basso prevede l’interazione. Modello senza interazione:

## 
## Call:
## occuMulti(detformulas = c("~CT.sens", "~CT.sens"), stateformulas = c("~Area+EL+SL", 
##     "~Area", "~0"), data = dt4, control = list(maxit = 1000, 
##     trace = TRUE, REPORT = 1), maxOrder = 2L)
## 
## Occupancy (logit-scale):
##                        Estimate    SE      z  P(>|z|)
## [ibex] (Intercept)       -0.682 0.307 -2.221 2.63e-02
## [ibex] AreaSB            -0.576 0.490 -1.175 2.40e-01
## [ibex] AreaSU             0.174 0.424  0.411 6.81e-01
## [ibex] AreaTB            -1.295 0.615 -2.106 3.52e-02
## [ibex] EL                 0.270 0.187  1.448 1.48e-01
## [ibex] SL                 0.368 0.169  2.172 2.98e-02
## [leopardo] (Intercept)   -0.753 0.305 -2.469 1.36e-02
## [leopardo] AreaSB         1.790 1.569  1.140 2.54e-01
## [leopardo] AreaSU         2.291 0.465  4.926 8.40e-07
## [leopardo] AreaTB        -2.142 0.789 -2.713 6.66e-03
## 
## Detection (logit-scale):
##                          Estimate     SE      z   P(>|z|)
## [ibex] (Intercept)         -3.520 0.1116 -31.55 1.93e-218
## [ibex] CT.sensmedium        0.684 0.2266   3.02  2.54e-03
## [leopardo] (Intercept)     -3.052 0.0661 -46.17  0.00e+00
## [leopardo] CT.sensmedium   -1.565 0.4402  -3.55  3.79e-04
## 
## AIC: 3924.488 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 109 
## Bootstrap iterations: 0

Modello con interazione:

## 
## Call:
## occuMulti(detformulas = c("~CT.sens", "~CT.sens"), stateformulas = c("~Area+EL+SL", 
##     "~Area", "~1"), data = dt4, control = list(maxit = 1000, 
##     trace = TRUE, REPORT = 1), maxOrder = 2L)
## 
## Occupancy (logit-scale):
##                             Estimate    SE     z  P(>|z|)
## [ibex] (Intercept)            -1.174 0.378 -3.11 1.87e-03
## [ibex] AreaSB                 -1.257 0.664 -1.89 5.84e-02
## [ibex] AreaSU                 -0.606 0.552 -1.10 2.73e-01
## [ibex] AreaTB                 -0.933 0.648 -1.44 1.50e-01
## [ibex] EL                      0.265 0.187  1.42 1.56e-01
## [ibex] SL                      0.374 0.169  2.21 2.71e-02
## [leopardo] (Intercept)        -1.407 0.427 -3.29 9.97e-04
## [leopardo] AreaSB              2.047 1.275  1.61 1.08e-01
## [leopardo] AreaSU              2.530 0.517  4.90 9.80e-07
## [leopardo] AreaTB             -1.876 0.814 -2.30 2.12e-02
## [ibex:leopardo] (Intercept)    1.484 0.560  2.65 8.01e-03
## 
## Detection (logit-scale):
##                          Estimate    SE      z   P(>|z|)
## [ibex] (Intercept)         -3.524 0.112 -31.47 2.53e-217
## [ibex] CT.sensmedium        0.686 0.227   3.02  2.51e-03
## [leopardo] (Intercept)     -3.051 0.066 -46.26  0.00e+00
## [leopardo] CT.sensmedium   -1.495 0.392  -3.81  1.39e-04
## 
## AIC: 3918.218 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 104 
## Bootstrap iterations: 0

4.5 Interazione leopardo - lupo

##         model negLogLike nPars      AIC   delta     AICwt
## 1 lup_leo_dip   1748.239    11 3518.478 0.00000 0.6472368
## 2 lup_leo_ind   1749.846    10 3519.692 1.21383 0.3527632

In questo caso il modello più supportato contiene l’interazione. Modello senza interazione:

## 
## Call:
## occuMulti(detformulas = c("~CT.sens", "~1"), stateformulas = c("~Area", 
##     "~EL+Dis", "~0"), data = dt5, control = list(maxit = 1000, 
##     trace = TRUE, REPORT = 1), maxOrder = 2L)
## 
## Occupancy (logit-scale):
##                        Estimate    SE     z  P(>|z|)
## [leopardo] (Intercept)   -0.725 0.303 -2.39 1.68e-02
## [leopardo] AreaSB         5.786   NaN   NaN      NaN
## [leopardo] AreaSU         2.263 0.464  4.88 1.07e-06
## [leopardo] AreaTB        -2.166 0.789 -2.75 6.05e-03
## [lupo] (Intercept)       -0.764 0.221 -3.46 5.35e-04
## [lupo] EL                -0.299 0.205 -1.46 1.44e-01
## [lupo] Dis                0.506 0.198  2.55 1.07e-02
## 
## Detection (logit-scale):
##                          Estimate     SE      z   P(>|z|)
## [leopardo] (Intercept)      -3.05 0.0661 -46.15  0.00e+00
## [leopardo] CT.sensmedium    -1.83 0.2355  -7.79  6.67e-15
## [lupo] (Intercept)          -3.95 0.1431 -27.61 8.13e-168
## 
## AIC: 3519.692 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 62 
## Bootstrap iterations: 0

Modello con interazione:

## 
## Call:
## occuMulti(detformulas = c("~CT.sens", "~1"), stateformulas = c("~Area", 
##     "~EL+Dis", "~1"), data = dt5, control = list(maxit = 1000, 
##     trace = TRUE, REPORT = 1), maxOrder = 2L)
## 
## Occupancy (logit-scale):
##                             Estimate    SE     z  P(>|z|)
## [leopardo] (Intercept)        -0.996 0.337 -2.95 3.16e-03
## [leopardo] AreaSB              1.375 1.020  1.35 1.78e-01
## [leopardo] AreaSU              2.265 0.465  4.87 1.10e-06
## [leopardo] AreaTB             -2.075 0.791 -2.62 8.68e-03
## [lupo] (Intercept)            -1.169 0.331 -3.53 4.15e-04
## [lupo] EL                     -0.231 0.205 -1.13 2.60e-01
## [lupo] Dis                     0.408 0.205  1.99 4.62e-02
## [leopardo:lupo] (Intercept)    0.738 0.450  1.64 1.01e-01
## 
## Detection (logit-scale):
##                          Estimate    SE      z   P(>|z|)
## [leopardo] (Intercept)      -3.05 0.066 -46.21  0.00e+00
## [leopardo] CT.sensmedium    -1.45 0.390  -3.73  1.95e-04
## [lupo] (Intercept)          -3.94 0.143 -27.62 6.00e-168
## 
## AIC: 3518.478 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 46 
## Bootstrap iterations: 0

4.6 Interazione stambecco - lupo

lup_ibex_dip <- occuMulti(detformulas = c('~1','~CT.sens'), stateformulas = c('~EL+Dis', '~Area+EL+SL', '~1'), data = dt1, control=list(maxit=1000, trace=TRUE, REPORT=1))
lup_ibex_ind <- occuMulti(detformulas = c('~1','~CT.sens'), stateformulas = c('~EL+Dis', '~Area+EL+SL', '~0'), data = dt1, control=list(maxit=1000, trace=TRUE, REPORT=1))
##          model negLogLike nPars      AIC     delta     AICwt
## 2 lup_ibex_ind   1230.308    12 2484.616 0.0000000 0.5418256
## 1 lup_ibex_dip   1229.476    13 2484.951 0.3353889 0.4581744

In questo caso il modello più supportato non contiene l’interazione. Modello senza interazione:

## 
## Call:
## occuMulti(detformulas = c("~1", "~CT.sens"), stateformulas = c("~EL+Dis", 
##     "~Area+EL+SL", "~0"), data = dt1, control = list(maxit = 1000, 
##     trace = TRUE, REPORT = 1), maxOrder = 2L)
## 
## Occupancy (logit-scale):
##                    Estimate    SE      z  P(>|z|)
## [lupo] (Intercept)   -0.764 0.221 -3.462 0.000537
## [lupo] EL            -0.299 0.205 -1.463 0.143355
## [lupo] Dis            0.506 0.198  2.552 0.010716
## [ibex] (Intercept)   -0.682 0.307 -2.221 0.026351
## [ibex] AreaSB        -0.576 0.490 -1.175 0.239802
## [ibex] AreaSU         0.174 0.424  0.411 0.681279
## [ibex] AreaTB        -1.295 0.615 -2.106 0.035181
## [ibex] EL             0.270 0.187  1.448 0.147610
## [ibex] SL             0.368 0.169  2.172 0.029825
## 
## Detection (logit-scale):
##                      Estimate    SE      z   P(>|z|)
## [lupo] (Intercept)     -3.950 0.143 -27.61 8.45e-168
## [ibex] (Intercept)     -3.520 0.112 -31.55 1.94e-218
## [ibex] CT.sensmedium    0.684 0.227   3.02  2.54e-03
## 
## AIC: 2484.616 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 107 
## Bootstrap iterations: 0
## 
## Call:
## occuMulti(detformulas = c("~1", "~CT.sens"), stateformulas = c("~EL+Dis", 
##     "~Area+EL+SL", "~1"), data = dt1, control = list(maxit = 1000, 
##     trace = TRUE, REPORT = 1), maxOrder = 2L)
## 
## Occupancy (logit-scale):
##                         Estimate    SE      z P(>|z|)
## [lupo] (Intercept)        -0.599 0.258 -2.323 0.02017
## [lupo] EL                 -0.290 0.208 -1.391 0.16419
## [lupo] Dis                 0.524 0.199  2.630 0.00855
## [ibex] (Intercept)        -0.520 0.331 -1.572 0.11587
## [ibex] AreaSB             -0.558 0.491 -1.135 0.25637
## [ibex] AreaSU              0.225 0.427  0.528 0.59766
## [ibex] AreaTB             -1.329 0.614 -2.165 0.03035
## [ibex] EL                  0.264 0.189  1.398 0.16209
## [ibex] SL                  0.367 0.169  2.169 0.03011
## [lupo:ibex] (Intercept)   -0.581 0.456 -1.274 0.20272
## 
## Detection (logit-scale):
##                      Estimate    SE      z   P(>|z|)
## [lupo] (Intercept)     -3.955 0.144 -27.55 4.37e-167
## [ibex] (Intercept)     -3.521 0.112 -31.53 3.46e-218
## [ibex] CT.sensmedium    0.683 0.227   3.01  2.61e-03
## 
## AIC: 2484.951 
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 98 
## Bootstrap iterations: 0

##Selezione dei modelli multispecie, a coppie, testando l’ipotesi di interazione

##           model negLogLike nPars      AIC      delta     AICwt
## 2  dom_ibex_ind   3319.739    17 6673.477 0.00000000 0.5093794
## 1  dom_ibex_dip   3318.776    18 6673.552 0.07504392 0.4906206
## 11  dom_leo_dip   3837.190    16 7706.380 0.00000000 0.7496418
## 21  dom_leo_ind   3839.286    15 7708.573 2.19340610 0.2503582
## 12  dom_lup_dip   3118.852    14 6265.704 0.00000000 0.7772587
## 22  dom_lup_ind   3121.102    13 6268.203 2.49952463 0.2227413
## 13   ib_leo_dip   1944.109    15 3918.218 0.00000000 0.9583079
## 23   ib_leo_ind   1948.244    14 3924.488 6.26971507 0.0416921
## 14  lup_leo_dip   1748.239    11 3518.478 0.00000000 0.6472368
## 24  lup_leo_ind   1749.846    10 3519.692 1.21383043 0.3527632
## 25 lup_ibex_ind   1230.308    12 2484.616 0.00000000 0.5418256
## 15 lup_ibex_dip   1229.476    13 2484.951 0.33538885 0.4581744

4.7 Riassunto grafico della selezione dei modelli per le interazioni inter-specifiche

Estimates for pairwise co-occurrence probabilities. Error bars show 95% confidence interval. Light blue indicates pairwise interactions supported by AIC, while orange indicates non supported pairwise interactions.

Estimates for pairwise co-occurrence probabilities. Error bars show 95% confidence interval. Light blue indicates pairwise interactions supported by AIC, while orange indicates non supported pairwise interactions.

Estimates for pairwise co-occurrence probabilities. Error bars show 95% confidence interval. Light blue indicates pairwise interactions supported by AIC, while orange indicates non supported pairwise interactions.

Estimates for pairwise co-occurrence probabilities. Error bars show 95% confidence interval. Light blue indicates pairwise interactions supported by AIC, while orange indicates non supported pairwise interactions.

4.7.1 predizioni multi-species

## $Livestock
## $Livestock$KS
##         mean standard_dev           SE 
##   0.60534397   0.14784248   0.07206244 
## 
## $Livestock$SB
##         mean standard_dev           SE 
##   0.44057281   0.15844494   0.08097191 
## 
## $Livestock$SU
##         mean standard_dev           SE 
##   0.63971750   0.25437865   0.07258252 
## 
## $Livestock$TB
##         mean standard_dev           SE 
##   0.29884588   0.12227894   0.07414067 
## 
## 
## $Ibex
## $Ibex$KS
##         mean standard_dev           SE 
##   0.34078750   0.10464736   0.08210319 
## 
## $Ibex$SB
##         mean standard_dev           SE 
##   0.20373393   0.06154862   0.06897559 
## 
## $Ibex$SU
##         mean standard_dev           SE 
##   0.37829737   0.11169770   0.08623092 
## 
## $Ibex$TB
##         mean standard_dev           SE 
##   0.15006472   0.06132802   0.06992785 
## 
## 
## $`Snow-leopard`
## $`Snow-leopard`$KS
##         mean standard_dev           SE 
##   0.32015771   0.00000000   0.06640912 
## 
## $`Snow-leopard`$SB
##         mean standard_dev           SE 
##    0.7382111    0.0000000    0.3053235 
## 
## $`Snow-leopard`$SU
##         mean standard_dev           SE 
##   0.82314430   0.00000000   0.05122325 
## 
## $`Snow-leopard`$TB
##         mean standard_dev           SE 
##   0.05237575   0.00000000   0.03617987 
## 
## 
## $Wolf
## $Wolf$KS
##         mean standard_dev           SE 
##   0.30501212   0.08089061   0.06422433 
## 
## $Wolf$SB
##         mean standard_dev           SE 
##    0.3444104    0.1217007    0.0682746 
## 
## $Wolf$SU
##         mean standard_dev           SE 
##   0.38605452   0.12440922   0.08421927 
## 
## $Wolf$TB
##         mean standard_dev           SE 
##   0.25220791   0.06903757   0.06117556